#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "LatoolException.h"
#include "IComplexMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class ComplexGeneralBandMatrix : public IComplexMatrix {
	private:
		long * ipiv;
		doublecomplex * ab;
		int m_rownum;
		int m_colnum;
		int m_ku;
		int m_kl;
		int m_ldab;
		bool m_isFactorized;

		ComplexGeneralBandMatrix(ComplexGeneralBandMatrix ^ orig) {
			m_rownum = orig->RowNum;
			m_colnum = orig->ColNum;
			m_ku = orig->m_ku;
			m_kl = orig->m_kl;
			m_ldab = 2 * m_kl + m_ku + 1;

			ab = new doublecomplex[m_ldab * m_rownum];

			for (int i = 0; i < m_ldab * m_rownum; i++) {
				ab[i] = orig->ab[i];
			}
			ipiv = new long[m_rownum];
			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}
			m_isFactorized = orig->IsFactorized;

		};
	public:
		ComplexGeneralBandMatrix(int rownum, int colnum, int subnum, int supernum) {
			m_rownum = rownum;
			m_colnum = colnum;
			m_ku = supernum;
			m_kl = subnum;
			m_ldab = 2 * m_kl + m_ku + 1;

			ab = new doublecomplex[m_ldab * m_rownum];

			for (int i = 0; i < m_ldab * m_rownum; i++) {
				ab[i].r = 0.0;
				ab[i].i = 0.0;
			}
			ipiv = new long[m_rownum];
			m_isFactorized = false;
		};

		~ComplexGeneralBandMatrix(){
			delete[] ab;
			delete[] ipiv;
		};

		virtual IComplexMatrix ^ Clone() {		
			return (IComplexMatrix^) gcnew ComplexGeneralBandMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property Complex default[int, int] {
			Complex get(int i, int j) {
				if ((i == j) || ((j > i) && (j - i <= m_ku)) || ((i > j) && (i - j <= m_kl))) {
					doublecomplex c = ab[j * m_ldab + (m_kl + m_ku + i - j)]; 
					return Complex(c.r, c.i);
				} else {
					return Complex(0.0, 0.0);
				}
			}
			void set(int i, int j, Complex value) {

				if ((i == j) || ((j > i) && (j - i <= m_ku)) || ((i > j) && (i - j <= m_kl))) {
					doublecomplex c;
					c.r = value.Real;
					c.i = value.Imag;
					ab[j * m_ldab + (m_kl + m_ku + i - j)] = c;
					m_isFactorized = false;
				} else {
					if (value.Real == 0.0 && value.Imag == 0.0) {
						return;
					} else {
						throw gcnew LatoolException("Value " + value.ToString() +  "doesn't fit to (" + i.ToString() + ", " + j.ToString() + ")");	
					}
				}
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::ComplexGeneralBand; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {
			throw gcnew LatoolException("Inversion for gengeral band matrix is not supported by lapack.");
		};

		virtual void Factorize() {

			long m = m_rownum;
			long n = m_colnum;
			long kl = m_kl;
			long ku = m_ku;
			long ldab = m_ldab;

			long info;

			zgbtrf_(&m, &n, &kl, &ku, ab, &ldab, ipiv, &info);

			m_isFactorized = true;
		};

		virtual void Solve(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char trans = 'N';
			long n = m_rownum;
			long kl = m_kl;
			long ku = m_ku;
			long nrhs = gb->ColNum;
			long ldab = m_ldab;
			long ldb = n;
			long info = 0;

			zgbtrs_(&trans, &n, &kl, &ku, &nrhs, ab, &ldab, ipiv, gb->dat, &ldb, &info);
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			long n = m_rownum;
			long kl = m_kl;
			long ku = m_ku;
			long nrhs = gb->ColNum;
			long ldab = m_ldab;
			long ldb = n;
			long info = 0;

			zgbsv_(&n, &kl, &ku, &nrhs, ab, &ldab, ipiv, gb->dat, &ldb, &info);

			m_isFactorized = true;

		};
	};
}